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The strongly correlated liquid state of a bilayer of charged particles has been studied via the HNC 
calculation of the two-body functions. We report the first time emergence of a series of structural 
phases, identified through the behavior of the two-body functions. 

PACS numbers: 71.45.-d, 73.61.-r, 64.70. Ja, 64.90.+b 

Electronic bilayer systems, consisting of two quasi-two-dimensional layers of electron (or hole) liquids, separated by 
a distance d comparable to the interparticle separation a (the Wigner-Seitz radius) within the layers, have revealed 
an unexpected richness of features and have been the object of many recent investigations. While most of the interest 
has originated from the quantum Hall effect, and thus the related studies have been directed at bilayers in strong 
magnetic fields, more recently it has been recognized that unmagnetized bilayers also constitute a physical system of 
remarkable complexity. Electronic bilayer systems in semiconductors can now be quite routinely fabricated by modern 
nano-technology in a fairly wide range of densities (i.e. r s values) and interlayer separations. Interestingly, virtually 
similar systems have come into existence also under quite different circumstances in ionic traps, where the charged 
particles are classical ions and the layers develop as a result of the system seeking its minimum energy configuration 
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At high enough r s {— a/ciB, as = is the effective Bohr radius) or T (— ) values (estimated to be between 
r = 98 and T = 222 Q ) the bilayer is expected to crystallize into a Wigner lattice. While there is no direct 
experimental evidence on the formation of such an electronic solid in semiconductor bilayers without magnetic field 
. (since the required r s values do not seem to be attainable by present-day experimental technique), the formation of 
£ — . ' layered crystalline structures has been observed in molecular dynamics computer simulations [^| and in ionic traps 
ON ■ |§. At the same time, a series of rather thorough theoretical investigations |j],|]-|6) have been carried out predicting 
the formation of classical Wigner lattices in electronic or ionic bilayers. These studies have revealed the existence of 
distinct structural phases, whose characteristics depend on the layer separation, i.e. on the d/a ratio: (see inset in 
Fig. 2): the staggered rectangular (ii); the staggered rhombic (iv), and the staggered triangular (v). Structure (ii) 
. comprises the simple triangular (i) and the staggered square (iii), as special cases. A similar quantum mechanical 
c"j ' study for the magnetized bilayer [|7| has led to very similar predictions. 

The understanding of the experimentally more important (and theoretically perhaps more challenging) liquid phase, 
on the other hand, is quite incomplete. Its analysis hinges upon the availability of the intralayer and interlayer two- 
body functions, (sn(r) and gi2 (r), respectively). While attempts have been made to determine the gij-s through 
various approximations, none of the results obtained can be considered reliable especially with regard to correctly 
reflecting the developing short range order: Zheng and McDonald || and Neilson and collaborators ||, have used the 
STLS (Singwi-Tosi-Land-Sjolander) approach, whose limitations are well-known; Kalman, Ren and Golden [To| have 
used a quasi-iterative method, which breaks down for small layer separations. 

The purpose of this Letter is to provide the first reliable calculation of the pair correlation functions gij (r) and 
related quantities such as the structure functions Sij(k) and the dielectric matrix £ij(k) for a bilayer system. The 
results show that different short range ordered structures develop already in the liquid phase, at relatively modest 
coupling values, emulating the different structural phases in the Wigner lattices. The gradual increase of the layer 
separation is accompanied by dramatic changes in g^ (r) , Sy (k) and (fc) that can be partially interpreted in terms of 
the structural changes referred to above, partially in terms of the new type of a substitutional orded-disorder transition. 
The bilayer seems to be a unique liquid system exhibiting structural transformations of the phase pursuant to the 
change of a system parameter: such transformations are commonly known only in the solid phase. 

Our method of approach is based on the classical HNC (Hyper-Netted Chain) approximation. This method has 
proven to be extremely reliable and accurate for Coulomb systems, both in three dimensions Jl2]Jl^] and two dimensions 
fl3| pL5j . We have adopted the HNC method to the bilayer situation, by mapping the (single component) bilayer onto 
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a two-component 2D system, with an interaction potential matrix [g_5 16 
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and using the two-component equivalent of the HNC method. The two layers are assumed to have equal densities n. 
The resulting system of equations for gijir) — hij(r) + 1, Cjj(r), and their Fourier transforms hij(k) and c%j{k) 

Mr) = e*M-<*«-*M - l ;M fc) = gf) ;Mt ) = ^W+^^W ; (2) 

[1 - ncnikJl 2 ~ [nc 12 {k)\ 2 l-ncn(fc) 

has been solved following Lado's original work on 2D systems [ fL5| . The underlying model is a classical one, where 
the average kinetic energy is represented by the inverse temperature parameter f3 and the intralayer coupling by the 
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parameter T = . This description is appropriate for bilayers arising in ionic traps if the confining potential is 
steep enough so that the degree of freedom perpendicular to the bilayer planes is suppressed; it is expected to be 
also quite adequate for the analysis of the spin-independent structural phases of the zero-temperature electron liquid, 
in the domain of high enough coupling (a > as) where the electrons are quasi-localized and the effect of intralayer 
exchange becomes small. Interlayer tunneling effects also disappear under the usually less stringent d > as condition. 
Under these circumstances, invoking the equivalence T ~ 2r s is a reasonable prescription for the interpretation of the 
results derived from the present model of semiconductor bilayers. 

Our result for the intralayer 311(f) and interlayer 312(7) two-body functions are summarized in Figs.l which show 
their variation as functions of the coupling T and of the interlayer separation d/a. In order to quantitatively analyze 
the information contained in the results we have used as principal indicators: 

(i) 312(0), the measure of the probability of finding a particle in layer 2 above another particle in layer 1: when the 
interlayer correlations become significant, 312(f) approaches zero; 

(ii) Ri, R11, Rin and Qj, Qn, Qui, the positions of the first three maxima of 311 (r) and 312(7"), respectively: 
these maxima identify the correlation-shells around a given particle, both in its own layer and in the layer adjacent 
to it; 

(iii) the intralayer coordination numbers pi, pn, pm and the interlayer coordination numbers ai, an, am, 
associated with each shell; the coordination numbers are defined by pi — 2rm J^' rgu(r)dr, ai = 2rm f^? rgi2(r)dr 

where M i (N { ) and M i (7V 4 ) are the minima preceding and following Ri{Qi). 

Short range order, signaled by the onset of oscillating behavior for 311(f) appears around T s = 3, more or less 
the same value as for the isolated 2D layer. The value of T s doesn't seem to be significantly affected by the layer 
separation d. Once the short range order has developed, the most dramatic novel feature of the two-body functions 
is the series of quite abrupt shifts in the positions and changes in the amplitudes of the first few peaks both of 311 (r) 
and of 312(f)- This can be well observed in Figs. 1 and 2. These features unambiguously point at the formation 
and phase-transformation-like changes of short range structures in the bilayer electron liquid, which can be paralleled 
with the results of recent studies of the bilayer crystalline phase fl],|]-§| ■ 

The physical features of the system are revealed through the two-body functions, the structure functions and the 
dielectric matrix of the system. In order to bring out the characteristic features in the details of the structure of 
two-body functions we have used the rather high T = 80 value for the latter. For the structure functions and for the 
dielectric matrix we have used the more realistic r — 20 and T = 10 values. 

1. Two-body functions 311(f), 312(f)- 

1.1 Interlayer correlations, as inferred from 311(0) seem to completely disappear for about d/a = 3, even for very 
high r values (see inset Fig. 2). Thus for layer separation higher than this value the two layers behave as virtually 
independent 2D systems. 

1.2 The positions, heights, and widths, and the resulting coordination numbers of the correlation-shells shift and 
change drastically as the interlayer separation is increased (see Fig. 2). These changes can be brought into cor- 
respondence with the expected formation of the different lattice structures in the solid phase as discussed above. 
Since each lattice structure carries a precise set of nearest neighbor positions {R a }, {Qa} and coordination number 
{p a }, {cr Q }(with a — 1,2,3... representing the first, second, third, etc. nearest neighbors in the lattice), these can be 
identified with the previously defined {R a }, {Q a }, {Pa}, and {a a } (values with a = I, II, III, representing the first, 
second and third shell positions). (The identification is not one-to-one: a given shell, in general, with the possible 
exception of the first shell, contains more than one set of nearest neighbors). Thus the formation of a series of "liquid 
phases" emulating the underlying lattice structural phases can be observed. 

1.3 At d = 0, when the two layers are collapsed into a single 2D system with double density 2n, their two two-body 
functions 311(f) and 312(f) are identical ,with Ru/Ri = Qn/Qi ~ 1.9 ratio, slightly higher than the characteristic 
•s/3 of a triangular lattice. The Ri — Qi — 1.25 value itself is somewhat below the lattice constant R\ = 1.34, 
similarly to what has been obtained in the Monte Carlo and HNC fl3|-|l5| calculations for an isolated 2D electron 
liquid. The double density triangular structure can also be regarded as a staggered rectangular lattice with a side 
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ratio 02/01 = The identity of the two two-body functions gn and g\2 and the accompanying equality of the 
coordination numbers, pi = <7j, however, clearly indicates that the two species are in a complete substitutional disorder, 
occupying the lattice sites at random. 

1.4 With d increasing away from the zero value, the substitutional disorder is rapidly replaced by a substitutional 
order appropriate for the staggered rectangular structure: by d/a — 0.5, pi and 07 assume their expected p\ = 2 and 
ai=4: values. At the same time, the rectangular unit cell deforms toward a square- like shape: this is well shown by the 
reduction of Ru/Ri to the ~ 1.5 value. The substitutional disorder that prevails in the domain < d/a < 0.5 would 
seem to require that vestiges of the peaks of g\\ (r) show up in #12(7) and vice versa. This, however, is apparently not 
the case: Ri and Qi, Rn and Qn, etc. are quite distinct. Nevertheless, the separations between Ri and Qi remain 
small enough and the peaks wide enough to accommodate the substitutionally ill-positioned particles. In fact, the 
Qui shell appears to be thriving on the ill-positioned population: by the establishment of the substitutional order 
the //-shell at Qn = 3.65 vanishes and is replaced by the more distant position at Qm = 4.99. 

1.5 The transformation to an underlying staggered square lattice structure becomes complete at d/a = 0.7. The 
indicators of this transformation are the jump of Rji from 2.30 to 3.89 and Rn from 3.87 to 5.63, and the accompanying 
jump of the coordination numbers pi from ~ 2 (same as p\ for the rectangular lattice) to ~ 7 (close to p\ + P2 = 8 
(pi = 4,p2 = 4) for the square lattice). All these changes reflect the higher symmetry of the square lattice. 

1.6 In the domain d/a > 0.7 the staggered square lattice gradually transforms into a staggered rhombic lattice with 
90 > (p > 60°: this allows the energetically favorable increase of the lattice constant, at the cost of the decrease of the 
energetically less important second and third neighbors. This trend is very clearly observable in the behavior of Ri, 
Ril, and Riu and also detectable in Qj, and Qui- There is no significant change in the coordination numbers in 
this domain, except a slight reduction in pi, pu, 07, an reflecting the reduced distance between the corresponding 
shells. 

1.7 Around d/a = 1.6, tp reaches the 60° value and in each layer the underlying structure becomes a standard 2D 
triangular (hexagonal) lattice, with lattice constant appropriate for density n. This can be inferred from Ri reaching 
its maximum and Rn, Rin and Qi, Q11, Qm reaching their minimum values, the former very close to \[2 times the 
Ri(d = 0), Ru(d = 0) and Rm(d = 0) values respectively. 

1.8 The triangular lattices emerging as the end product of the rhombic structure are staggered in such a way that 
the vertices of layer 2 lie over the midpoints of the sides of the triangles in layer 1 (and vice versa). This is not the 
minimum energy configuration for well-separated layers: this latter is attained when the vertices of the 2-triangles lie 
over the centers of the 1-triangles. Transition to this final configuration, consisting of a rigid translation of lattice 2 
with respect to lattice 1, takes place in the domain 1.5 < d/a ~ 2.0 and is shown by a slight, but perceptible increase 
in Qi and Qm. 

2. Structure functions Sii(fc), Si2(k) Fig. 3). The expected small k behavior of the structure function is governed by 
a generalized compressibility sum rule [jl7|: at k = the Sii(0) = — 512(0) condition is required; both of the values of 
5ii(0) and of dS g 1 ^ (k = 0) are governed by the quantity Lu —L12, the difference between the inverse compressibility 
and the inverse trans-compressibility: Lij = ^-K a , where K a , P and A are respectively the compressibility of the non- 
interacting gas, the pressure and the surface area. At d = the purely correlational L12 compensates the correlational 
contribution to the Lu (resulting in 5ii(0) = —512(0) = 1/2). This is also the consequence of the perfect screening 
sum rule which requires that n J d 2 r{hu(r) + h%2(r)} = —1. For high d values Lu becomes dominant, approaching 
the value appropriate for the isolated 2D layer, which, for T > 1 is Lu = 1 — 0.821T. The inset in Fig. 2 shows the 
L\\ — L\i values, extracted from Su(k) and S\2(k), corroborating the expected behavior. 

3. Dielectric response matrix. The inverse of the static dielectric response matrix rj(k) = e~ 1 (k) can be expressed 
in terms of the structure functions as T}ij{k) = 5ij — f3nipu(k)Sij(k). It determines the total potential <!>i in the two 
layers when either of them is perturbed by an external potential <f>i : $i = rjij $ j . 

At small layer separation the effective dielectric function deteij has the features of the dielectric function of an 
isolated 2D layer (see Fig. 4) with the characteristic anti-screening behavior for k < fc* ( ~ 25). This behavior 
becomes qualitatively different precisely at the layer separation where Lu — L12 changes from positive to negative 
(see inset in Fig. 4) Jl7| : en and £12 develop alternatively anti-screening and screening domains. 

We can compare our results concerning the critical values of d/a at the phase transition boundaries with those 
of Goldoni and Peeters (GP) [Quito our d/a values Table I shows the comparison. The d/a [(ii) — ► (Hi)] value has 
been taken along the liquid-solid phase boundary at the lowest V value (r = 118) in Fig. 8 of their paper. GP claim 
a wide range of stability (0.79 < d/a < 1.10) for the staggered square lattices. We don't see this happening: the 
gradual transformation into a rhombic structure seems to follow immediately the formation of the square lattice. The 
expected first order character of the (iv) - (v) transition is masked in the present description. Neither is there a clear 
indication of the transition taking place at f> = 69.6° predicted by GP) rather than at <p = 60° the indicator would 
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be a slight decrease in the pi values accompanying the detectable slight increase in the Qi values). The possibility 
for substitutional disorder is not part of the GP lattice model which allows for structural changes only, although, as 
our analysis reveals, it would be, for a finite temperature, a determining factor of the phase structure of the system 
in the d/a < 0.5 domain. 

In summary, we have calculated, the two-body functions, structure functions and the static dielectric matrix of a 
strongly correlated electronic bilayer liquid through the classical HNC approach. We have found a series of dramatic 
changes in the liquid structure as the layer separation varies from d/a = to d/a = 3. These changes can be brought 
into correspondence with the structural transformations and with substitutional order-disorder transformations in the 
underlying lattice structure. This behavior seems to be unparalleled in other strongly correlated liquid systems. We 
have also shown how the structure functions and the static dielectric matrix reflect these structural changes. 
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visiting at the University of California at San Diego and is grateful to Tom O'Neil for his hospitality. We would like 
to thank Dan Dubin for the discussions and for making his unpublished data available and to Ken Golden of the 
University of Vermont for a series of useful conversations. 
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CAPTIONS FOR FIGURES. 

Fig. 1. Variation of the structure of the two-body functions with increasing layer separations d/a at T = 80: (a) 
9u{r); (b) 0i2 (r). 

Fig. 2. Characteristic variation with layer separation of the first, second and third maxima of <?n(r) and g 12(f)' 
(a) positions Ri, Ru(b) associated coordination numbers pn, pu, pm and 07, an, 07//. The insets show the five 
principle lattice structures and the value 512 (^) at r = 0. 

Fig. 3. Variation of the structure functions with increasing layer separations d/a for r = 20 for Su(k) and Si2(fc)- 
Note that the coordinates of Sn and Si 2 are shifted with respect to each other. 

Fig. 4. Diagonal and off-diagonal elements of the static dielectric matrix for different layer separations (a) d/a = 
0.05, (in > L12) (b) d/a = 0.1 (in < L12) at T = 10. The insets show deteij and in — L12, the difference between 
diagonal and transinverse compressibilities. 

CAPTION FOR TABLE 
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Phase boundaries ((i): simple triangular; (ii/a): substitutionally disordered staggered quadrangular; (ii/b): substi- 
tutionally ordered staggered quadrangular; (iii): staggered square; (iv): staggered rhombic; (v): staggered triangular) 



Table: 



Phases 




d/a 






This work 




GP 


(ii/a) — ► 


- 0.5 






(ii) — > (in) 


0.7-0.8 




0.79 


(w) -> (u) 


1.5 - 2.0 




1.87 
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